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ABSTRACT 



Context. Modelling the emission properties of compact high energy sources such as X-ray binaries, AGN or 7-ray bursts 
represents a complex problem. Contributions of numerous processes participate non linearly to produce the observed 
spectra: particle-particle, particle-photon and particle- wave interactions. Numerical simulations have been widely used 
to address the key properties of the high energy plasmas present in these sources. 

Aims. We present a code designed to investigate these questions. It includes most of the relevant processes required to 
simulate the emission of high energy sources. 

Methods. This code solves the time-dependent kinetic equations for homogeneous, isotropic distributions of photons, 
electrons, and positrons. We do not assume that the distribution has any particular shape. We consider the effects of 
synchrotron self-absorbed radiation, Compton scattering, pair production/annihilation, e-e and e-p Coulomb collisions, 
e-p bremsstrahlung radiation and some prescriptions for additional particle heating and acceleration. 
Results. We illustrate the code's computational capabilities by presenting comparisons with earlier works and some 
examples. Previous results are reproduced qualitatively but some differences are often found in the details of the 
particle distribution. As a first application of the code, we investigate acceleration by second order Fermi-like processes 
and find that the energy threshold for acceleration has a crucial influence on the particle distribution and the emitted 
spectrum. 
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1. Introduction 

High energy sources, such as X-ray binaries, active galactic 
nuclei (AGN hereafter), or 7-ray bursts, exhibit spectra de- 
tectable to very high energy. This radiation must originate 
in a plasma for which a significant fraction of the particles 
have relativistic energies. Understanding the properties of 
these hot plasmas remains a challenge in the modelling of 
X- and 7-ray sources. 

Among the many processes at work, there are particle- 
particle interactions, such as Coulomb collisions, particle- 
photon interactions such as Compton scattering, syn- 
chrotron radiation, bremsstrahlung emission, or pair pro- 
duction/annihilation, and particle- wave interactions that 
lead to particle acceleration. However, the way they add or 
compete is highly non-linear and the cross sections involved 
are complex. Investigating a large parameter space is re- 
quired, and in spite of important breakthroughs, these plas- 
mas are still poorly understood. Analytical studies provided 
interesting qualitative results with approximations, but a 
more general approach based on numerical simulations is 
required to explain the details and complexity of contempo- 
rary observations. The first detailed investigations were an- 
alytical attempts to model the Co mpton scattering in ther- 
mal plasmas of fixed temperature (iBisnovatvi-Kogan et all 
ll97lUSunvaev fc Titarchuklll980l : lGuilbertlll98lHZdziarskil 
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119851 : ICuilberti [t986). In parallel, some of these re- 
sults were confirmed by Monte-Carlo simulations (e.g. 
iPozdnvakov et all 119831: iGorecki fc Wilczewski|[l98l . The 
additional role of pair production and annihilation in ther- 
mal plasmas, whose temperatures were determ ined self- 
consistently, was then studied both anal y tically (ISvenssonl 
Il982bl 1198.4 iGuilbertfcStepnevl fl985l : iKusunosd Il987h 
and numerically (|Zdziarskil 119841 , Il985f ) . These works con- 
stituted significant advances because they explicitly ac- 
counted for the back reaction of the radiation field on 
the plasma temperature. However, they were limited to 
thermal distributions of particles, whereas significant ev- 
idence of strongly non-thermal populations was found in 
many sources. For instance, spectra of blazars or radio 
loud AGN were shown to be shaped at least by the syn- 
chrotron s elf-compton emission of purely non-thermal elec- 
trons fe.g. iGhisellini et al. 19983). At these high energies, 
accelerated particles cool on very short timescales before 
they can be thermalized for instance by two body colli- 
sions. The balance between this cooling and acceleration 
typically produces non-thermal distributions. Acceleration 
processes are still poorly understood. A simple way to sim- 
ulate the effect of particle acceleration is to inject parti- 
cles at high energy. Although it does not reproduce ex- 
actly the physics involved, this prescription has been widely 
used and produced interesting results (as shown in most 
of the references cited here). Significant effort has also 
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been taken in developing more precise modelling of ac- 
celeration mechanisms, but in such studies, the radiation 
field is treated crudely ([Li et al.lll996l: iDermer et al.lll996t 
Hi fc Millerlll997t iKatarzvhski et al.ll2006bh . 

With the increasing number of considered processes 
and the increasing precision of their description, numeri- 
cal analysis has become a prime method of investigation. 
Even so, a full treatment of the problem accounting for the 
coupled evolution of inhomogeneous, anisotropic distribu- 
tions of leptons and photons, both in momentum and po- 
sition spaces, appears to be still beyond the capabilities of 
present-day computers. Numerical simulations of high en- 
ergy plasmas have been performed mainly following three 
different approaches which all make trade-offs between the 
various aspects of the problem. 

Fi rst, the Monte C arlo technique (jPozdnvakov et aTl 

1 19771 1 Stern et alj 1995) allows one to follow particles and 
photons in space, time, and energy as they undergo mu- 
tual interactions. It solves the full radiative transfer prob- 
lem and accounts explicitly for geometrical effects. At 
present, the MC method is probably the mos effective 
way to model fully 3-dimensional problems. However, this 
detailed procedure is time consuming, particularly when 
modelling the rapid dynamics o f the non-thermal electro n 
population in momentum space (jMalzac fc Jourda in 2000) , 
and when synchro tron self-absorptio n effects are important 
(see discussion in [Stern et al.l [i~995). For this reason, the 
Monte-Carlo methods have been applied to date to pure 
Maxwellian plasmas and/or steady state problems with 3D 
geometry. 

Another method that accounts correctly for the geome- 
try, involves solving numerically the exact radiation trans- 
fer equation for given geome tries and particle distributions 
(jPoutanen fc SvenssonlfT996l ). This method is far more ef- 
ficient than Monte Carlo simulations which makes it eas- 
ier to compare with data. It is, however, far less versatile 
than Monte Carlo methods and does not solve the kinetic 
equations for particles. The back reaction of the radiation 
field on the particle distribution is modelled only for the as- 
sumption of a Maxwellian plasma (in which case the plasma 
temperature may be adjusted according to energy balance). 
The method applicability is also limited to the resolution 
of steady state problems. 

The third approach, which we adopt in this paper, aban- 
dons the detailed description of the geometry to concen- 
trate on the kinetic effects. It consists of solving the lo- 
cal kinetic equations for the particle and photon distribu- 
tions. To maximazc efficiency, radiative transfer is usually 
modelled with a simple photon escape probability formal- 
ism assuming isotropic photon and particle distributions. 
This method can be applied to different, possibly time- 
dependent, problems for which geometry does not play a 
crucial roleQ- Within the limits of the one-zone approxima- 
tion, it is more efficient than other methods and allows for 
fast data fitting. 

The first detailed investigations of high energy plasmas 
with this technique concentrated on thermal pair plasmas 
(|Fabian et all Il986t iGhisellinil [l987h . More precise mod- 



1 For problems in which geometry is important, radiative 
transfer can in principle be accounted for by coupling this ki- 
netic code with a radiati on transfer solver o r a M onte Carlo 
code (as demonstrated by iBottcher k, Liand (|200lf )h although 
computing time may then become a serious issue. 



elling was then proposed in which the particle distributions 
were decomposed into the sum of a t hermal low-energy pool 
and a n arbitrary high energy tail (Lightman & Zdziarski 
1 19871: ISvenssonl 1987t ICoppil 119921 [Zdziarski et all 119931 : 
iGhisellini et"afl 119931: |Li_eLalJ 1199J). The latter models 
have been applied most to fitting and interpreting data. 
They do not however describe the possible deviation from 
a Maxwellian distribution at low energy, nor do they ad- 
dress explicitly any thermalization process. Only recent 
numerica l work considered fully arbitrary distributions of 
particles. IGhisellini et al.l (|l998al) concentrated on the role 
of synchrotron self-absorbed radi ation in AGN. They c on- 
firmed previous analytical results (jGhisellini et al.lll988f ) by 
demonstrating that the exchange of energy between parti- 
cles by means of synchrotron photons can be an efficient 
thermalization process in magnetized sources. These sim- 
ulations focused however on this specific interaction, and 
other processes were only considered by crude approxima- 
tions, particularly Comptqn rad iation, or not considered 
at all. iNavakshin fc MeTial (|l998f) investigated the thermal- 
ization of arbitrary distributions by two-body particle in- 
teractions and heating by high energy protons. The addi- 
tional role of synchrotron radiation was not however con- 
sidered. The most complete numerical treatment of high 
energy plasmas was p robably one developed i n the con- 
text of 7-ray bursts bv lPe'er fc Waxmanl (|2005t ). Our code 
is similar to this study but differs in that these authors 
considered neither particle stochastic acceleration nor the 
effect of Coulomb losses 

The code presented here solves the time-dependent 
equations simultaneously for isotropic, arbitrary photon, 
electron, and positron distributions. The evolution of these 
populations is modelled in time while being affected by 
self-absorbed cyclo-synchrotron radiation, Compton scat- 
tering, pair production/annihilation, e-e and e-p Coulomb 
collisions, self-absorbed e-p bremsstrahlung radiation, and 
additional particle acceleration and heating. Each process 
is described with minimal approximations and by using in 
most of cases the exact cross sections. For instance, the 
formulae used for the synchrotron emission and absorption 
are valid from the sub-relativistic to the ultra-relativistic 
regime. This numerical strategy allows one to investigate 
many different astrophysical situations that occur in vari- 
ous high energy sources. 

The structure of this paper is as follows. Section [2] pro- 
vides a description of the microphysics adapted into our 
code. Then, in Sec. [3J we present the numerical techniques. 
Finally, in Sec. [3] the code is tested against previous pub- 
lished results, providing an overview of its capabilities. 

2. Radiation and kinetic processes 

We describe the processes included in the code. We present 
first the general notation used in this paper. The parti- 
cle energy is described cither by the relativistic Lorentz 
factor 7 = E/m e c 2 , by the adimensional momentum p = 

V/m e c = y / 7 2 — 1, or by the beta parameter [5 — p/7, 
where m e is the electron mass and c is the speed of light. 
Similarly, the photon energies are described by their fre- 
quencies v or oj — hv/m e c 2 . The particle and photon pop- 
ulations are described by their angle- integrated distribution 
functions iV e ± = and N v = Jn^;, where dAf e ± and 

dN v are the number of electrons, positrons and photons per 
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unit volume d 3 x and per unit momentum dp or frequency 
dv. For simplicity, the total lepton distribution is also used: 
N e — N e - + N e +. Finally, R is the typical length scale 
of the emission region. Since we consider a homogeneous 
medium, most of the processes are scale-free, meaning that 
most quantities are simply proportional to R, R 2 , or R 3 . 
For those quantities, R determines only the overall normal- 
ization factor. For instance, the total luminosity of unmag- 
netized sources scales as i? 3 , but there is no reference scale 
in the problem. The only process that explicitly involves 
a reference length scale is the synchrotron self-absorption, 
since it is independently determined by both the magnetic 
intensity and the total magnetic energy of the source, which 
is related directly to the source size for a given magnetic 
intensity. 

2.1. Self-absorbed Synchrotron Radiation 

The cyclo-synchrotron radiation process is produced by 
charged particles gyrating along magnetic field lines. It is 
one of the most important processes in astrophysics and is 
invoked to explain the radio emission of many magnetized 
sources, such as supernova remnants, pulsar wind nebulae, 
AGN, X-ray binaries, and 7-ray bursts. In particular, it pro- 
duces soft photons that can be up-scattered by Compton 
scattering producing the well known double humped syn- 
chrotron self-compton spectra used to model for example 
the emission of blazars or radio-loud AGN. The reverse 
process, cyclo-synchrotron absorption, is less well known, 
although efficient at low energy. Besides explaining some 
observed spectra, cyclo-synchrotron emission and absorp- 
tion both influence the involved particles by cooling and 
heating them respectively. As we discuss later, these inter- 
actions can therm alize high energy part icles (the so-called 
synchrotron boiler, IGhisellini et alJll988f l. 

Following the main assumption of the code, we assume 
isotropy in the radiation field and particle pitch angle. 
These are accurate approximations when the magnetic field 
is tangled. The cyclo-synchrotron emission and absorption 
are characterized by the emissivity spectrum j s (p, v) (erg 
s _1 Hz _1 ) of one single particle of momentum p and the 
cross section a s (p,v) (cm 2 ) Both quantities are related to 
each other by the formulgn (iLe R oux 1961; Ghiscllini et al.l 
1988; IGhisellini fc Svenssonlll991h : 



1 



4-7T 2m e v 2 p 2 



dp \pjj s (p, v)} 



(1) 



where m e is the electron mass. The emissivity and cross sec- 
tion depend on the magnetic field, whose i ntensity is char- 
acter ized by the magnetic compactness (|Ghisellini et alJ 
EUl): 



I, 



a T R B 2 
m e c 2 8tt 



As mentioned earlier, the synchrotron self-absorption de- 
pends explicitly on the source size. Although the overall 
normalization of the emissivity and absorption is only pro- 
portional to a combination of the magnetic field intensity 
and source size (namely Ib), their shape depends on the cy- 
clotron frequency vb = eB/2irm e c, which depends only on 



2 The term 1/4-7T in Eq.[T]and|8]relates to the average over the 
solid angle 



the magnetic field intensity. For a given magnetic compact- 
ness parameter, simulations of sources with different sizes 
correspond to cases with different magnetic field intensities 
and therefore produce different observed spectra. 

In uniform systems, the time evolution in the mean in- 
tensity integrated over solid angles I v — hvcN v (erg s^ 1 
Hz -1 cm~ 2 ) is described by the equation: 



with 

pb v = / N e j s (p,v)dp 



(erg cm 3 s 1 Hz ) , 



k v = / N e a s (p 1 v)dp (cm 1 ) . 



(3) 



(4) 
(5) 



We note that this equation differs from one often used i n 
previous works (for example IGhisellini et al]|l988l Il998a! ). 
These have concentrated mostly on the steady state prop- 
erties of magnetized sources and did not include this time 
dependence. To account for photon escape and the non- 
absorbed part of the observed spectra, a finite size do- 
main is assumed in these papers and the space equation 
d x Iv = (J*v — k u I u is solved for uniform emissivity fj, s and 
absorption k s on a typical length scale R, which yieds the 
synchrotron self-absorbed radiation: I u — /i/k(1 — e~ KR ). 
This is an approximate way to deal with the space de- 
pendence of the simulated system since it implies a non- 
uniform synchrotron radiation whose feedback onto the 
lepton equation would need to be incorporated in a fully 
space-dependent model. This solution would holds also 
only when synchrotron self-absorption is involved. Other 
processes, such as Compton scattering or pair produc- 
tion/annihilation, would also contribute to the emissivity p, 
and absorption n in some way and in this description it is 
unclear how the synchrotron interaction couples with others 
in the global equation for the radiation field. We consider 
a homogeneous (or averaged) radiation field which, asso- 
ciated with the photon escape probability (see discussion 
in Sect. 12. 6|) represents another approximate way to model 
crudely the geometry. However, this method solves the ex- 
act time-dependent equation and the synchrotron emission 
is added consistently to other emission processes. 

Simultaneously, the equ ation for the t i me evolution of 
the l epton populations is ijMcCravl 119691 : IGhisellini et alJ 
1988): 



d t N r 



AlN n 



D B .N n 



(6) 



where the j/p factors result from the choice of p as a vari- 
able instead of 7 and 



(2) At = 



1 



Dl = 



mc* 
1 1 



(j S -0- s Iy) dv 

3 ah 



47r m 2 c 2 



-dv 



(7) 
(8) 



Equation [6] can be written in many different ways that are 
analytically equivalent. The one used (Eq. [6][8]), associated 
with a specific numerical scheme to estimate derivatives, 
enable good numerical accuracy. In particular, the energy 
conservation can be easily satisfied to machine precision 
when particles and the photon field exchange energy, since 
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dE e /dt = mc 2 J 7(9 P [A e N e ] dp = - JJ {j s - <r s I v )N e dpdv = 
— J{p s — n s Iv)dv = —dE v jdt. 

There is no exact analytical expression for pitch angle 
averaged- and photon direction integrated- emissivity and 
absorption that is valid in all regimes. The exact values 
are derived from numerical integrations, which are time- 
consuming and hard to perform, especially for low energy 
particles when the emission is dominated by a few narrow 
harmonics. However, approximations have been proposed, 
which are valid in some regimes (e.g. iMarcowith fc Malzad 
120031) . We use a combination of two approximations. For 
sub-relativisti c particles, we u se the formula for j s first 
proposed by IGhisellini et alj (|1998af) and corrected by 



Kat arzvhski et al 



(|2006aF both to match the relativistic 
spectrum more accurately and to describe the spectrum 
more closely close to the minimal frequency. This approx- 
imation is less accurate for very low energy particles (typ- 
ically [3 < 0.1). However, for most astrophysical cases, the 
emission is produced mainly by energetic particles, and this 
regime has little influence on the total particle distribution 
and radiation field. Numerical experiments have confirmed 
that the choice of j s and a s at low energy has negligible 
effect. In the relativistic regime, we use the well-known 
synchrotron power spectrum integrated over an isotropic 
distribution of pitch angles (jCrusius fc Schlickeiserl Il986t 
IGhisellini et al.|[l988h . We note that we correct this for- 
mula similarly to match more accurately the sub-relativistic 
regime. The transition between both regimes is achieved by 
applying an exponential threshold/cut-off to the formulae 
at 7 = 2. The formulae for o~ s are computed analytically 
from j s with Eq. [1] and implemented into the code. 



2.2. Compton Scattering 

Compton scattering is a well-known interaction between 
leptons and photons. Most of the previous studies assume 
thermal distributions of particles. In such cases, only the 
temperature and total number of particles are computed. 
This assumption enables rapid computation with simple 
formulae but omits the physics of non-thermal particles. We 
therefore use exact analytical expressions for unpolarized 
radiation and arbitrary distributions of isotropic particles 
and photons. 



2.2.1. Basic equations 

The effect of Compton scattering can be described by the 
sum of individual encounters over the entire distributions. 
The scattering of isotropic photons of energy hvo = u>om e c 2 
off isotropic particles of energy Eq = jomc is fully char- 
acterized by the resulting distribution of scattered photons 
o'cCpoj^'o ~~ * v )- This differential Klcin-Nishi na cross sec- 
tion has been computed by several authors (|Jonesl [T968t 
lBrinkmannlll984tlNagirner fc Poutanenlll994D . The numer- 
ical evaluation of these analytical expressions can be dif- 
ficult, in particular for low or high energy particles and 
photons. We use an expression based on the formulae by 
iJoned (I1968D and mod ified to overcome numerical accuracy 
issues (|Belmontll2008h . 

The exact time evolutions of the full particle and photon 
distributions are described by the following equations: 



9tN e ±(p) 



-N e ±{p) J N„{is )ca c (vo,p)dv , (9) 

d t N v {v) = JJ ' N e {p )N u (v )ccr c {p 0l VQ -> v)dp dv 

-N v {v) J N e (p )c<j c Q {v,p )dpo , (10) 

where the photon frequency v{p) is constrained by the en- 
ergy conservation during one scattering event: hv{p) — hvo + 
jm e c 2 — 7om e c 2 = 0. For each distribution, the first integral 
provides the number density of scattered particles/photons 
that have a particular energy after one single scattering and 
the second one is the probability that particles/photons of 
this energy are scattered to some other energy. This is what 
we refer to as the integral approach. As discussed hereafter, 
the numerical computation of this integral suffers from ac- 
curacy issues because of discretization. 

In the small-angle scattering limit, when the scattered 
photons (or particles) have energies similar to those incom- 
ing, a Fokker-Planck approximation can be used (FP here- 
after). In this case, a second-order series expansion of the 
exact equations produces the FP evolution equations for 
the different species: 



d t N v {v) = [A C V N V ] + ^ [D C V N V ] 



d t N e ±( P ) = d, 



, t -AlN e± 



ld p -D c N e ± 
p \p 



(11) 
(12) 



with 



A c , 



Al = 



N e ca\(jp,v)dp , Dl 



N e ca c 2 (p,v)dp , (13) 



N v cal(p,v)dv , D c e = / N v c(jl{p,v)dv , (14) 



and where we have introduced the first 3 moments of the 
scattered photon distribution: 



CTo(Po,^o) 



c c (po, — > v)dv , 



(15) 



ctJ(po,^o) = / (w -u )o- c (po,v — > v)dv , (16) 
^(Pa, v Q ) = / (w - ujo) 2 o- c (po, v n -> v)dv , (17) 



N e ±(p )N u (is )co- c {po,is Q -> v{p))dp Q dv 



which are the total cross section, the mean photon energy, 
and the dispersion, respectively. 

This approximation allows far quicker computation 
since, when the first moments have been tabulated, only 
single integrals are required whereas the exact computa- 
tion requires double integrals. However, the Fokker-Planck 
approach used to model the evolution of particle and pho- 
ton distributions is valid only in regions of the incident 
energy space (votPo) for which the relative energy ex- 
change in one scattering is small: Ap(po,vo)/po << 1 and 
Av(pq, vo)/v{yo) « 1, respectively. Thes e conditions wi ll 
be presented in a forthcoming publication (|BelmonlJl2008l ). 

2.2.2. Numerical strategy 

In contrast to the Fokker-Planck approximation, the inte- 
gral approach is exact analytically. However, when used to 
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compute the evolution numerically, it leads to some numer- 
ical issues directly relate d to the use of non-linear grids 
(jNavakshin fe MelialH99l . 

With logarithmic grids, the energy bin size is larger at 
high energy. For example when, low energy photons are up- 
scattered from high energy particles, their relative energy 
gain is high, and these photons are scattered numerically 
from low energy bins to higher energy bins. During this in- 
teraction, the particles lose only a small fraction of their 
energy. If the energy bin size is too large, these particles 
remain in their original bin and numerically, they do not 
lose energy. The energy balance is not therefore exactly sat- 
isfied and the error can propagate and become large when 
the density of low energy photons is also high. Although 
less relevant to most astrophysical situations, a symmetrical 
problem appears when high energy photons are scattered 
by low energy particles. 

This numerical issue is not present in regions of the 
incident energy space (pq,vq) for which the scattered 
distributions are far wider than the energy bin size: 
Ai/(p , vq)/8v{v ) » 1 and Ap(p , is )/pa » 1 for the 
evolution of the photon and particle distributions respec- 
tively. After selecting the ranges and resolution of the pho- 
ton and particle energy grids, these conditions constrain the 
region in which the integral approach is valid. Fortunately, 
the regions for which the integral and the Fokker-Planck 
approaches are valid are in part complementary. The code 
therefore combines the two approaches: 

• For the particle evolution: 

In the integral approach, the integration over the photon 
distribution in Eq.[5]is only completed above a given photon 
energy v c {po) that depends on the incident particle energy, 
whereas the integrals on frequency in Eq. [H] are completed 
up to v c in the Fokker-Planck approach. The total time 
evolution is then given by the sum of both contributions: 
d t N e± = (d t N e± ) Fp + (d t N e± ) lntcgial . 

• For the photon evolution: 

A similar combination is used for the photon equation, with 
the definition of a critical particle momentum p c (i>o) under 
which the FP approach is used and above which the integral 
approach is used. 

If the number of energy bins per decade is too small, 
the validity domains for both calculations may become in- 
dependent and the accuracy of the computation may de- 
crease. This will be true, however, only for small regions 
of the grids for which there are few particles and photons, 
corresponding to very small errors. By combining the two 
approaches, we find empirically that an energy resolution 
of typically 10 energy bins per decade provides errors that 
are smaller than other truncation errors. 



2.3. Pair production/annihilaton 

As for Compton scattering, we describe the pair produc- 
tion and annihilation for the case of isotropic distributions 
of particles and photons. Single photon-photon pair pro- 
duction and pair annihilation events are characterized by 
the differential cross sections <r p (y\ , v 2 — > p) which corre- 
sponds to the pair (i.e. electrons or positrons) momentum 
spectrum produced by the recombination of photons of fre- 
quencies v\ and v 2 , and <7 a (p_,p+ — * v\ which corresponds 
to the emission spectrum generated by the annihilation of 
one electron of momentum p_ and one positron of momen- 



tum p + , respectively. Then, the evolution of the distribu- 
tions is described by: 



N u (v 1 )N u (v 2 )c(J v {v 1 , v 2 ^p)dv 1 dv 2 
-N e± {p) J N^(p')a%(p,p')dp' , (18) 

N e -(p-)N e +(p + )ccr a (p-,p + -*■ v)dp^dp + 
-N v {v) [ N„{v')a%(v,v')dv' (19) 



d t N e ±(p) 



d t N„{v) 



where, as for Compton scattering, the zeroth mo- 
ment of both the annihilation spectrum CTq(p_,]3 + ) = 
1/2 J cr a (p_,p+ — > v)dv and the pair-produced distribu- 
tion <Jq{vx,V2) = 2/ o v (v\tV2 — * p)dp have been usec0. 
The analytical expressions for photon-photon pair produc- 
tion and pair annihilation cor respond to Eqs. (24-29) of 
Boettcher fc Schlickeiserl (|199<D and Eqs. (23,33,55-58) of 
SvenssonT i 1982al ). respectively. 

In contrast to Compton scattering, there are no numer- 
ical problems in computing directly the integral over the 
particle and photon distributions, even for low resolution 
grids. 



2.4. Coulomb scattering 

Two kinds of Coulomb-type interactions are considered: 
scattering of leptons by other leptons and scattering of 
leptons by protons. When the particles are not too ener- 
getic, e-e collisions tend to thermalize the pair distribu- 
tions. In some astrophysical situations, protons are assumed 
to have a high temperature, so that e-p collisions tend to 
heat the lepton populations. Both kinds of interactions are 
described by the Boltzmann collision integral. Computing 
this integral numerically is challenging, mainly because the 
Coulomb cross section diverges when the energy exchange 
becomes too small. We assume instead the approximation 
of small angle scattering, which leads to simple Fokker- 
Planck equations. As has been already discussed in the liter- 
ature, the contribution of large angle scattering is negligible 
(permer fc Liandll989t INavakshin fc Melialll998fh 



2.4.1. Moeller and Bhabha e-e scattering 



The FP coefficients (A e ~ e and D e 



for Moeller e~ 



and Bhabha — e scattering are similar. The relative 
difference is typically of the order of ~ 1/lnA, where th e 
Coulomb logarithm In A is large (|Dermer fc Liand 119891) . 
Neglecting these terms, the FP coefficients are computed 
from the following integrations over the mirror distribu- 
tions: 



K± e iP±) = J N eT(p T )a e (P-,P+)dp T , 
D e e± e {p±) = J N eT (p T )d e (p-,p+)dp T . 



(20) 
(21) 



The specific coeffic i ents a e and d e were first given by 
INavakshin fc Melial ([1998) (Eq. 24 and 35 respectively). 



3 The 1/2 and 2 factors result from the fact that one pair anni- 
hilation produces 2 photons and one photon-photon annihilation 
produces 2 leptons. 
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The equation for d how e ver co ntains typos that were cor- 
rected by Eq. 6 of HE! (|2003). 



2.4.2. Coulomb e-p scattering 

In some astrophysical situations, protons are believed to 
have temperatures far higher than those of electrons. In 
these cases, Coulomb collisions with leptons tend to heat 
the electrons. As for the e-e collisions, the FP coefficients 
(A e ~ p and D e ~ p ) for e-p Coulomb scattering are computed 
from integrals over the proton distribution: 



N p (p p )a p (p e ,pp)dpp 



D e± (Pe) = / N p (p p )d p {p e ,Pp)dpp 



(22) 
(23) 



where a D and d v are derived from Eqs. (45-48) of 
iNavakshin fc Melial (|1998l ). The code could calculate the 
exact proton distribution as for electrons and positrons. 
However, it would add a fourth kinetic equation and re- 
quire more computational time. We use instead a thermal 
proton distribution. Depending on the physical situation 
being modelled, the proton temperature can be set to equal 
a constant at the beginning of the simulations or evolved 
with time to provide a constant electron heating (see sec- 
tion 



2.5. Particle and photon injection 

The code also allows for injection of particles into the sys- 
tem. This can represent a real injection from an outer source 
(e.g. particles from the standard disc into an ADAF). But 
injection of high energy particles is most commonly used 
to mimic particle acceleration processes. Any distribution 
N™t can be injected at each time step. Thermal, Gaussian, 
power-law, and mono-energetic injections have already been 
included into the code. The injection of particles is con- 
trolled by the particle injection compactnesifl: 



L 



R 2 <T T Air 



(24) 



Similarly, the code allows for photons injection that can 
account for example for seed photons from the cold accre- 
tion disc in X-ray binaries. The photon injection rate is 
controlled by the parameter: 



R 2 a T 47r 



>R/a T 



hv 

m e c 2 



N^dv 



(25) 



So far the code can account for photon injection with a pure 
black-body spectrum of specific temperature or a multi- 
temperature black-body spectrum characterised by the in- 
ner and outer temperatures. 

2. 6. Particle and photon losses: geometry of the source 

In general, photons and particles can also escape from the 
system. The precise way in which they escape depends on 



4 Note that some authors use instead the kinetic energy to 
define the compactness parameter: l e ± = 4-kR 2 (tt /3c J(-y — 
1)TO P . 



the detailed geometry of the simulated source, which goes 
beyond the scope of such a 1-zone kinetic code. Although 
the losses must occur at the boundaries of the simulated 
plasma, we use a standard method and we consider all 
photons (or particles) to have the same averaged probabil- 
ity p° sc (or p e *±) of escape. We assume spherical geometry 
and use probability laws describing this geometry approxi- 
mately. The total luminosity of the source at each frequency 
is then: 



(26) 



The code allows for fully trapped pairs in extremely 
magnetized systems (no loss) and for freely escaping pairs 
(with an escape probability proportional to their velocity: 
P^ s ± = /?c/ R). Other escape laws can be defined and imple- 
mented easily. 

The photon escape is more debatable. The photon 
dynamics arc affected strongly by Compton scattering, 
high energy photons do not scatter and can escape freely, 
whereas, when the optical depth is large, low energy pho- 
tons can be scattered so significantly that they become 
trapped in the system. Depending on their energy, the ex- 
act way in which they escape strongly involves geometri- 
cal effects. In the code, we use the escape rate r^ sc (or es- 
cape probability defin e d as jp f„ sc = r° sc x R/c) derived by 
iLightman fc Zdziarskil (|198 



l/T v 

I CSC 



c/R 



l + rH/H/3 
where T* sc is the averaged escape time, 



t(uj — hv/m e c 2 ) = Rut 



N f-^2A d p 

(JT 



(27) 



(28) 



is the Compton interaction probability of photons (of fre- 
quency v) with the lepton distributions and 

r 1 for lj < 0.1 

f(ui = hv/m e c 2 ) = 1(1- w)/0.9 for 0.1 < w < 1 (29) 







for u> > 1 



is a relativistic factor correcting for the fact that forward 
collisions are less efficient in trapping the photons. 

The choice of escape probability is important and dif- 
ferent laws can lead to substantially different result^. 
Although it was shown that this escape probability repro- 
duces well the results of Monte Carlo simulations in a spher- 
ical g eometry (jLightman et al.lll 987; Lightm an fc Zdziarskil 
1987), the conclusions of IStern et al.l (|1995l ) imply that the 



escape probability may be slightly underestimated. Since 
the escape luminosity must equal the injected power in a 
steady state, a smaller escape probability implies a stronger 
radiation field inside the source. The consequences on the 
shape of the photon and lepton distributions become sig- 
nificant only when pair production and annihilation are ex- 
tremely efficient (typically at high optical depths). Figure 
[1] shows spectra in these cases when the deviation from the 



5 For instance, when com p aring h is res ults with those of 
ILightman fc Zdziarskil i|1987l ). ICoppil (|1992f ) attributed the dif- 
ference to different descriptions of the microphysics, whereas, 
from the results of several simulations, we believe that the dif- 
ference is due to a different choice for the escape rate: he used 
(J?/c)/(l + r(w)/(w)) instead of Eq.l27l 
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Monte Carlo simulations become significant. Other escape 
probabilities were proposed (e.g. IStern et "aT1ll995h to re- 
produce more successfully the results from MC simulations 
in some specific regimes, but none were shown to be fully 
consistent with MC results. The use of an escape probabil- 
ity to mimic geometrical effects is of course the main lim- 
itation of our code. However, significant deviations appear 
only in optically thick plasmas when steep gradients in tem- 
perature and intensity appear, whereas jets and coronae in 
XRB and AGN are optically thin media, the largest optical 
depths observed being r ss 2 — 3. The precise geometry of 
the sources is also unknown and describing accurately the 
escape probability in one peculiar geometry is therefore not 
necessarily helpful. 

2. 7. Additional particle heating/acceleration 

Particle heating and acceleration are probably amongst 
of the most mysterious problems of high energy sources. 
Observations show evidence for hot plasmas or high energy 
tails in the particle distributions, but little is known about 
the precise mechanisms that generate these populations. 
Most previous work did not address this problem directly. 
Non-thermal high energy particles were instead injected 
into the system with an arbitrary (usually power law) dis- 
tribution. This ad hoc injection assumes an instantaneous 
acceleration of particles. It does not take into account the 
fact that particle acceleration has to compete with other 
cooling processes. Another simple approach, often used to 
account for lepton heating, consists of assuming that power 
is provided by some unspecified process to the supposedly 
thermal distribution of electrons. These prescriptions for 
particle acceleration and heating arc implemented in the 
code. However, in addition, we also attempted to follow 
a more physical approach by implementing two additional 
specific mechanisms for heating and acceleration, namely 
Coulomb heating and second order Fermi acceleration. 

• e-p Coulomb-like heating: 
As has already been discussed, collisions with hot protons 
can heat the pair distributions. The way in which the inter- 
action is adapted into the code is described in Sect. 12.4.21 
When the Thomson optical depth is lower than unity, this 
heating is known to become inefficient and other processes 
must operate, which are not fully understood. A possible 
means of accounting for this additional heating is to mimic 
the he ating by thermal protons but with enhanced effi- 
ciency (jNavakshin fc M clia 1998). Although we do not aim 
to describe any precise microphysics, this heating prescrip- 
tion estimates consistently both FP coefficients: the heat- 
ing rate and its related diffusion coefficient. For this heating 
prescription, the temperature is set and the total number 
of protons is constrained by the initial neutrality. The usual 
cooling and diffusion coefficients A^± p (Eq. |22|) and D^± p 
(Eq. I23[) are then multiplied by an efficiency factor r\. This 
efficiency is computed at each time step so that the total 
heating is controlled by a constant heating parameter 



lc 



in R 2 ot 



- v A e I p N e dp 



(30) 



• 2nd order Fermi-like acceleration: 

This type of acceleration could be generated for example by 
the interaction between the electron and wave turbulence. 



Diffusion of particles in momentum space is then described 
by the general equation: 



df 
dt 



l_d_ 

p 2 dp 



p 2 D 



difr 



df 
dp 



(31) 



where f(p) is the phase-space density. When considering 
an equation about N e ± it yields a Fokker-Planck equation 
with the two coefficients: 



PI \ 7 

7 



(32) 



(33) 



We assume a Fermi-like process for particles with an en- 
ergy above some minimal energy, and we use D dlS — 
p2 e -~{pc/p) /2t acc , where p c is the threshold momentum, a is 
the threshold width (we use typically a = 3), and t acc is the 
typica l acceleration time of the particles (jKatarzvnski et aTl 
2006b). The FP coefficients are then: 



5 + 4p 



2 / Pc 

a 7 — 

P 



-(Pc/P)° 



= _L^L-(p<=/pr 



t -v 2 



,(34) 



(35) 



The precise values of acceleration time and the threshold 
frequency depend on the microphysics and turbulent prop- 
erties of the plasma, which are poorly known. We set in- 
stead the total energy injected into accelerated particles by 
defining a constant compactness parameter: 

4-7T R 2 <7 T 



-A^N e dp 



(36) 



and compute the corresponding acceleration time at each 
time step. 

2.8. Bremsstrahlung emission 

The bremsstrahlung process has several contributions to 
the system evolution: it produces additional soft photons 
that can then be up-scattered by high energy particles, it 
cools down emitting, high energy particles and, in the ab- 
sorbed part, it heats low energy particles. In arbitrary plas- 
mas, there are three different contributions: lepton-proton 
(e —p), electron-electron or positron-positron (e — e), and 



electron-positron (e 



") bremsstrahlung. 



Electron-proton self-absorbed bremsstrahlung is in- 
cluded in the code. Proton are assumed to have non- 
relativistic temperature and to be at rest in the plasma 
frame and the emission is generated by the motion of lep- 
tons in the external electrostatic potential of protons. The 
situation is formally the same as in the case of synchrotron 
emission, which is generated by the motion of leptons in 
an ex ternal magnetic field, and a similar formalis m can be 
used (|Le Rouxlll960HGhisellini fc Svenssonlll99lh . The ex- 
act interactio n cross section a® u v alid for all lepton energies 
(|Heitlerlll95i Uauch fc Rohrlichlll976l ) is used to compute 
the emissivity spectrum j ep (erg s" 1 Hz -1 ) of individual 
electrons and the absorption cross section a cp (cm 2 ): 

j ep (p,u) = hv{3col v {p, V )Nl ot (37) 
a ep (p,v) = —^—^—d p \pyj ep (p,u)] (38) 
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where N*" 1 is the proton density, (3 = v/c is the lepton 
velocity, and a ep is evaluated numerically. Approximations 
of the total loss rate and the spectrum emitted by a ther- 
mal plasma are recovered by integrating j ep over a thermal 
distribution of leptons. Then, the evolution equations for 
particles and photons are derived exactly as in the case of 
synchrotron emission (Eqs. • 

Electron-proton bremsstrahlung is the dominant con- 
tribution in low energy, e-p plasmas. In low-energy, pair 
plasmas the e + — e~ process dominates, whereas at high 
temperatures the major contribution originates in e — e 
bremsstrahlung. Differential cross sections for e — e and 
er — e + bremsstrahlung can be found in various regimes, 
either in the rest frame of one of the leptons or of the cen- 
tre of mass of the two interacting particles (jHeitlerl Il954t 
lAlexanianl 1196 8; Haugj [19751 Il985[ ). However, there is no 
formula in the frame of the plasma for the cross section 
integrated over all directions of the emitted photon; sim- 
pler thermal approximations are also often irrelevant since 
the cross section typically increases with particle energy 
and non thermal emission of high energy particles often 
dominates the overall bremsstrahlung emission. In princi- 
ple, these difficulties could be overcome numerically. 

However, in many astrophysical cases bremsstrahlung 
emission is insignificant. For plasmas in thermal equilib- 
rium, simple approximations were proposed for the cool- 
ing rates, which allowed for comparison with ot her pro- 
cesses (e.g. lGouldlll980t IStepnev fc Guilbertlll983f) . It was 
found that bremsstrahlung cooling dominates over pair- 
annihilation cooling and Coulomb relaxati on only of highly 
relativistic te mperatures: /csT > 1 MeV (|Svensso n 1982b; 
iStepnevll 19831 respectively). By integrating the synchrotron 
and Compton cooling rates over a hot thermal distri- 
bution of particles {UbT = 1 MeV) of radiation energy 
density U v , we find that, for plasmas with optical depth 
of the order of unity, bremsstrahlung emission dominates 
only when R<jTU v /m e c 2 + Ib < 4 x 10~ 3 , that is for 
unmagnetized, photon-starved plasmas. For these com- 
pactness parameters, only hotter plasmas have a signifi- 
cant bremsstrahlung contribution, although, for astrophys- 
ical sources, these high temperatures are unrealistic since 
pair production and annihilation tend to prevent tempera - 
tures reaching above a few hundreds keV (|Svensson| [l984). 
Similarly, it is has been shown that non-thermal particles 
emit primari ly synchrotron radiation, even for weak m ag- 
netic fields (|Wardzihski fc Zdziarskill200ri ICoppilll992f ). 

For these reasons, e — e and e~ — e + bremsstrahlung have 
not been included in our code and modelling of unmagne- 
tized, highly relativistic plasmas with a weak radiation field 
is postponed to future work. 



3. Numerical methods 

Including all processes outlined above, the total physi- 
cal system is described by the following set of 3 integro- 
differential equations: 



d t N v 
d t N e ± 



s,,± 



1 



+ [A V N V ] + -c£ 2 [DM > (39) 



P e ±N e ±+d p fr/pA e ±N e ±] 



ld p (lD e ±N e ± 
P \P 



(40) 



The source terms S v (t,v,N v ,N e - ,N e +) and 
S e ±(t,p, N v , N e -, N e +) combine the contributions of 
injection, Compton scattering (treated in the integral 
approach), annihilation/production, bremsstrahlung, and 
synchrotron emission: 

S v = N^{t,v) 

p OO pOO 

+ dv I dp N e (p )N l ,(v )ca c (po,vo;^) 

JO Jpa("o) 

+ // N e -(p^)N e +{p + )ca a (p-,p + - 1 v)dp-dp + 



(41) 



+ / dp dv N e (p )N v (v )ca c (po,vo;v(p)) 
Jo JMpo) 

+ JJ^ N v (y 1 )N v {v 2 )ca p (v 1 ,v 2 ;p)dv 1 dv 2 . (42) 

The loss terms Pv(t, v, N v , N e - , iV e +) and 
P e ±(t,p, N v , N e ^) also combine contributions from es- 
cape, Compton scattering, pair production/annihilation, 
Bremstrahlung, and synchrotron absorption: 

/•OO 

P„ = pT (v,N e -,N e+ ) + N e (p)ca c (p, v)dp 

J Va (v) 



S e ± = N^(t,p) 



;± J ( 

oo c OO 



+ / N v (v')ca%(v,i/)dv' + J N e c(a s + a cp )dp , (43) 

pOO 

P e ± =Plf(p)+ / N v (v)co%(p,v)dv 

Jv a (pa) 

+ ( N eT (p')ca a (p,p')dp' . (44) 



The total Fokker-Planck coefficients are then the sums of 
the individual coefficients defined for each process: 



A v = - \ N e (p)cal(v,p)dp 
Jo 

D v = I N e {p)ca 2 {v,p)dp 
Jo 



(45) 
(46) 



and 



A„± = 



1 



((is + iep) - c(a a + <7 cp )hvN v ) dv 



m e c 

+ / N u {v)ccr c 1 {v,p)dv 
Jo 



+ I N eT (p')a e (p,p')dp' + I N p (p')a p (p,p')dp' 



+ A c e ± p + AIT , (47) 



D„± = 



1 h f { j s + j op )N„ 



Air mic 



v 



r"c(p) 

+ I N v (v)ar%(v,p)dv 
o 



+ J N e +(p)d e (p,p)dp + j N p (p)d p (p,p)dp 

+ D^ p + D^ c . (48) 
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In this section, we describe the numerical strategy used to 
solve these equations. 



3.1. Tables 

Solving Eqs. [39] and [40] for the aforemen- 
tioned physical processes involves signifi- 
cant application of the many cross sections 
(is j ®s , Jcp j Ccp j o~ c j °~o i ) C2 , cr a , Cq , CTp , Oq , a e , d e , dp, dp). 
The Compton differential cross section cr c (z/ ,Po — ► i/) is 
a three-entry table that contains typically over tens of 
millions of elements, as is also the case for the differential 
cross sections for pair production/annihilation. These 
coefficients are evaluated once at the beginning of the 
simulation and stored into tables. This enables faster 
computation, although the memory requirements can 
become significant as the grid resolution increases. For 
instance, resolutions of above 256 points per grid require 
more than 100 Mo of RAM to store only one of these 
tables, which can be a limiting factor for some desktop 
computers. 



3.2. Boundaries 

The total set of equations to be solved is given by an 
integro-differential system. To account for particles and 
photons that have energies beyond the limits of the grids 
(very low- or very high energy particles/photons), specific 
conditions must be set at the grid boundaries uj min , u) max , 
p m in, and Pmax- For the differential Fokker-Planck part of 
the equation which corresponds to a local operator, the 
boundary conditions set values only for the ghost bins just 
behind the boundaries and are used to define the derivative 
at the boundaries. For the integral part, they include the 
physics of all particles and photons outside the grids. We 
have chosen to use wall- type boundary conditions. These 
conditions do not allow particles/photons to travel in and 
out of the grids and conserve their total number precisely. 
In the Fokker-Planck part, it corresponds to a zero- flux 
condition. For the integral part, a specific derivation of the 
differential cross section is completed: the total probabil- 
ity that particles/photons are scattered or produced out- 
side the grids is summed and added to the probability that 
they are scattered or produced at the final bin of the rel- 
evant boundary. As a result, all particles/photons remain 
inside the grids. These conditions are of course artificial 
and can generate spurious effects at the boundaries. While 
most particles/photons have energies inside the grids, these 
effects will, however, remain small. 

For example, although the total number of parti- 
cles/photons is conserved, the no flux condition for the 
FP part of the equations introduces a small energy flux 
at the boundaries which can be evaluated easily, that is 
the energy losses are dE u /dt = m e c 2 [DuiVu]""" and 
dE e ±/dt = m e c 2 [j/pD e ±N e ±]p^ for photons and lep- 
tons respectively. When the grid is not sufficiently large, 
this effect can introduce significant errors. By selecting the 
boundaries far from the bulk of the distribution however, 
the distributions N u and N e ± vanish and the losses are 
negligible. 



3.3. Numerical solver 

Solving the time-dependent problem is challenging for a 
number of reasons. 

First, the problem involves many different scales of en- 
ergy and time spanning many orders of magnitude. This 
implies subtracting very large numbers or multiplying very 
small numbers by very large ones, which can lead to nu- 
merical accuracy issues. 

This also involves considering very short timescales. For 
instance, when low energy photons are scattered by high 
energy particles, they absorb a significant amount of en- 
ergy instantaneously, which must be modelled with very 
small time steps. If the problem was linear and differential, 
the maximal time step required to guarantee the conver- 
gence and stability of an explicit scheme would be set b y 
the Courant condition (|Courant. Friedrichs fc Lewy| [T928). 
In that particular case, our simple scheme for the photon 
equation would be stable if St < min {Slu/A„; 2(Suj) 2 /D^} , 
where Slo is the bin size and the minimum is computed over 
the entire grid. Similarly, a condition depending on the mo- 
mentum bins would be set for the stability of the equation 
for particles. When logarithmic grids are used, the time step 
is set to a small value by the small bins at the low energy 
part of the grids. The equation is far more complicated, and 
there is no mathematical justification for using the Courant 
condition. However, the main idea remains that when the 
grids decline to low energy, the time step required to make 
an explicit scheme stable, quickly becomes too small to fol- 
low the evolution on the dynamical timescale R/c. 

In these cases, implicit schemes would be more effi- 
cient, since they are always stable. Implicit schemes are 
easy to implement and efficient only for local, linear prob- 
lems. When the problem is linear, the solver only has to 
inverse a matrix. For local problems such as differential 
ones, the matrix is sparse and rapidly inverted. However, 
the problem is highly non-local. For example, the integral 
approach of Compton scattering describes events in which 
photons in some energy bin can be scattered to some dis- 
tant bin by a single interaction. As a result, the evolution 
in the photon distribution at some energy is governed not 
only by neighbouring bins, but by the entire grid. The cor- 
responding matrix is dense and its inversion becomes more 
time consuming. In addition, the problem is highly non- 
linear, so that there is no exact implicit solution to such a 
problem. 

Considering these previous remarks, we implemented a 
semi-implicit scheme. By keeping the other distributions 
fixed, we solved the equation for the distribution N of one 
species (photons, electrons, or positrons) in the following 
way: 

JV"+ 1 - N n , , 

_ _ gn _ pn jy-n+1 

St 

+d x [A n N n+1 ] + ifl£ a [D n N n+1 ] (49) 

where the exponents n and n + 1 indicate the distributions 
and coefficients at two consecutive time steps. The scheme 
is not fully implicit but the equations for iV™ +1 , N™^ 1 and 
N™^ 1 can be solved easily by inverting a simple multi- 
diagonal matrix. This scheme only conserves the number 
of particles/photons and their energy to truncations errors. 
As explained earlier, using grids spanning many orders of 
magnitude introduces accuracy issues, which become severe 
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when number and energy are only conserved to truncation 
errors. For that reason, the code iterates the 3 equations 
represented by Eq. [49] at each time step until the conserva- 
tion laws are satisfied to some specific precision. 

Two spatial schemes were impleme nted to solve Eq. 
49l t he upwind Chang-Cooper method ([Chang k, Cooperl 
1970) and a more straightforward method that we devel- 
oped. The former is only accurate to first order in space, 
but was shown to contain the most superior numerical 
properties for solving Fokker-Pla nck equations with a few 
choices of constant coefficients ([Park fc Petrosianl [T996) ■ 
Compared with higher order schemes, it is more diffusive, 
i.e. more stable but less accurate. When radiative pro- 
cesses are involved, accuracy is important, and therefore 
we decided to use a second-order accurate scheme to 
estimate the derivative^]. This scheme is based on the 
use of two energy grids for each distribution (namely the 
centres and the faces of the bins) and the derivatives are 
computed as follows: {d x f) i = (f i+ i/ 2 - /i-i/a)/^ and 

(plf)i = ((/i+l - fi)/dXi+l/2 ~ ifi ~ fi-i)/dXi_ 1/2 )/dxi. 
All quantities at the bin boundaries are computed by 
linear interpolation: fi+1/2 = (fi + /i+i)/2, except the 
ip/l)i+i/2 factor in Eq.[T2]for particles, which is computed 
separately as (jpi + Pj+i)/(7j + 7;+i) to ensure accurate 
energy conservation. This simple centred scheme conserves 
the total number of particles and photons, and the total 
energy more precisely, although is less stable. All examples 
presented were performed using the latter scheme. 



3.4. Computation time 

Tests and applications shown in this work were com- 
pleted with medium energy resolution for which n v = 
128 — 256 and n p = 256, i.e. more than 20 bins per 
decade in particle momentum and more than 10 bins per 
decade in photon frequency. Compared w ith former codes 
(|Lightman fc Zdziarski|[l987HCoppilll992D that assume low 
energy particles are purely thermal, our code typically 
solves the equations for twice as many bins of particle mo- 
mentum. In addition, as presented, it adopts very few ap- 
proximations and only when they are valid. The code is 
therefore slower than some previous codes. All results pre- 
sented were derived with desktop computers with 1GHz 
processors and 512 MB RAM. Calculations have duration 
typically of between a few seconds and one hour, the most 
time-consuming step being the computation of the multiple 
integrals of the Compton and pair production/annihilation 
cross sections over the distributions. 



4. Tests and applications 

We present a first few applications of the code to check 
its capabilities and illustrate the problems that can be ad- 
dressed. Most examples presented consist of comparisons 
with previous work, although we also address a few more 
issues, such as the threshold for particle acceleration. 



The scheme is accurate to second order only when the grid 
is linear. Using a logarithmic grid as we do here actually reduces 
the accuracy. Nevertheless, numerical experiments with the code 
have shown that this scheme is more accurate than the Chang- 
Cooper method. 



4.1. Model with external soft photons 

Significant effort has already been expended on studying 
steady state solutions of unmagnetized sources. We only 
reiterate some known results in checking the code compu- 
tational capacitie s. As a first case, we reproduced the re- 
sults of Fig. 1 in ICoppil (|1992| ) with parameters typical of 
AGN. This case was modelled by injecting mono-energetic 
e + — e~ pairs at 7 = 10 3 and soft photons as a black body 
of temperature fc^T = 1.07 x 10~ 5 rn e c 2 . In this unmag- 
netized source, the length scale is unimportant. All effects 
were included except additional particle acceleration and 
e-p Coulomb scattering and e-p bremsstrahlung radiation. 
Figure ([T]) shows spectra obtained with the code and com- 
parisons with previous work. Examples of output are also 




i(r 6 io _ * io~ z 10° io ! 

cj = hi//m t c 2 



Fig. 1. Photon spectra for models with external soft pho- 
tons (l e = l e - + l e + = 1,10,100,1000 from lower to 
higher curve and l v — 2.5Z e ). Solid lines: this work; dashed 
lines: results with the EQPAIR code; and dotted lines: 
results of Monte Carlo simulation for spherical geometry 
([Stern et alj [l995). The spectra from EQPAIR were nor- 
malized to match the other ones. To simplify the c ompar- 
i son, the escape probability p c J c = 1/(1 + rf) of ICoppil 
(| 1992T ) is used in this figure. 



presented in Table [T] The temperatu re presented i n this ta- 
ble was computed using Eq. (2.8) of lCoppil l| 1992ft only for 
the thermal part of the particle distribution. The results 
are fully consistent with those computed with the latest 
version of the public EQPAIR code. The major deviations 
appear at the annihilation line for large optical depths and 
luminosities (l e = 100 and 1000). 

When T e /(l e + l v ) > 1, effects of pair annihilation and 
Coulomb scattering become significant. We investigated 
this regime by r eproducing the results presented in Fig. 
2 of lCoppil (|l992l ). This simulation had the same input pa- 
rameters as previously, except that particles were injected 
with a power law distribution (7 m i n = 1.4, 7 max = 10 3 , 
r = 2.4, l e — 1) and l v = 0.03. The particle distribution 
and photon spectrum for such a case are plotted on Fig. [2] 
The corresponding output parameters are listed in Table [2] 
Again, good agreement is achieved, particularly for the par- 
ticle distribution, confirmi ng that, in t his peculiar case, the 
approximations made bv lCoppil (|1992[ ) were valid. The only 
difference occurs at high energy in the photon spectrum. 
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L 


work 


PY» 


T . 
1 e 






0^2 — 10 




Coppi92 


1.7 


0.047 


7.9 




0.637 


1 


EQPAIR 




0.078 


12.6 


0582 


0.614 




This work 


1.70 


0.059 


8.64 


0.0611 


0.633 




Coppi92 


23 


0.502 


5.7 




0.863 


10 


EQPAIR 




0.508 


8.27 


2.70 


0.854 




This work 


22.6 


0.548 


4.50 


3.00 


0.901 




Coppi92 


87 


3.34 


2.2 




0.979 


10 2 


EQPAIR 




3.17 


2.56 


63.9 


0.996 




This work 


80 


3.22 


1.96 


61.9 


1.02 




Coppi92 


120 


12.4 


0.62 




1.42 


10 3 


EQPAIR 




11.95 


0.567 


644 


1.39 




This work 


112 


12.01 


0.44 


628 


1.41 



Table 1. Models with external soft photons: comparison 
with pre vious results. PY 3 is 10 3 times the pair yield de- 
fined in ICoppil (|!992f l. r e = Ra T J(N e - + N e +)dp is the 
total Thomson optical depth. 8 3 = 10 3 x kBT/m e c 2 is the 
temperature of the thermal part of the distribution. lx /l e 
is the ratio of the X-ray luminosity in the 2-10 keV band to 
the injection compactness parameter. ot2— to is the spectral 
index in the same energy band. 





PY 3 




#3 


lx/le 


CK2-10 


Coppi92 


3.54 


0.644 


306 


0.077* 


0.726 


This work 


3.27 


0.624 


273 


0.0733 


0.717 



Table 2. Output parameters for l e = 1 and l v — 0.03 . Same 
as ta ble [T] *Note that the luminosity obtained by ICoppil 
(1992) was multiplied by 10 to correct what we think is a 
typo. 



4.2. e-p Coulomb like heating 

We investigated the effect of coulom b-type heating and 
compa red the results with those of iNavakshin fc Melial 
(|l998l ). We consider an unmagnetized source heated by a 
Coulomb-like process. We assumed a closed system (no in- 
jection nor loss of particles) with a black-body soft-photon 
injection and studied its evolution under the effects of 
Compton scattering, e-e Coulomb exchange, pair produc- 
tion/annihilation, and e-p Coulomb-like heating. The pro- 
ton temperature characterising the final process was set to 
be 20 MeV. Two different cases were considered, the pa- 
rameters of which are given in Table [31 Regardless of the 







L 


r°- 


e„ 


Model 1 


420 


420 


0.05 


10~ 4 


Model 2 


8.4 


2.1 


0.02 


3 x 10~ 5 



Table 3. Input parameters for the runs on the study of 
the e-p Coulomb like heating. r°_ is the initial electron 
Thomson optical depth. Q„ = ksT /m e c 2 is the tempera- 
ture of the black body used for soft photon injection. 



initial particle distribution, it always evolved to the same 
steady solution. The thermalisation time however depended 
upon the precise initial s hape. Distributions in the t ransient 
phase were presented by INavakshin k, Melial (I1998D and we 




derived similar results. Figure [3] shows the steady distribu- 



P - Pi 

Fig. 2. Photons spectrum (upper panel) and particle distri- 
bution (lower panel) for power-law particle injection (with 
h = 1, 7min = 1.4, 7 max = 10 3 , T = 2.4) and for l v = 0.03. 
The resul ts with the c ode are shown in solid line and the 
results of ICoppil (|1992D are in da shed lines. Sin ce the par- 
ticle distribution at low energy in ICoppil (|1992D is assumed 
to be thermal but not resolved, we extended it down to the 
grid boundary with a Maxwell distribution for comparison. 



tions and spectra for both cases, and compares the result 
with previous work. The steady states correspond to quasi- 
thermal distributions. Both the c-c Coulomb exchange and 
the e-p Coulomb-like interactions thermalize the plasma ef- 
ficiently. In addition, the e-p Coulomb- like interactions heat 
the particle distribution. 

Our steady spectra are fully consistent with previous 
results: the Compton orders have the same amplitude, and 
the high energy spectrum breaks at the same energy. Some 
output parameters are also given in Table [H Our results 
confir m qualitatively the results by INavakshin k Melial 
(1998). The steady distributions have properties similar to 
those of a Maxwell-Boltzmann distribution but are nar- 
rower. Our results however correspond to systematically 
hotter distributions and a larger optical depth in the most 
energetic case (model 2). 

The heating efficiency coefficient is quite large: r\ « 
10 4 — 10 5 . This emphasizes the inefficiency of the e-p 
Coulomb collisions for example in cases typical of Seyfert 
galaxies. When the Thomson optical depth is far smaller 
than unity, they should be several orders of magnitude more 
efficient to reach the required heating rate. The efficiency 
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0.1 1.0 10.0 

P = P7 



Fig. 3. steady state spectra (upper panel) and electron 
distribution (lower panel) for e-p Coulomb-like heating. 
Solid and dashed lines sh ow the results of this work and 
iNavakshin fe MeTial (|l998f ). respectively. Input and output 
parameters are listed in Tables [3] and [3J respectively. 



4.3. Models with synchrotron soft photons 

Magnetized models have an additional source of soft pho- 
tons, which is the synchrotron emission from high energy 
particles. We investigated this case by studying generic 
cases of magnetised sources with no external source of soft 
photons and compared our results with those presented in 
Fig. 4 of ICoppil <| 19921 ) . Particles were injected at high en- 
ergy (7 = 10 3 ). We assumed they were trapped by the 
magnetic field and did not escape. The equilibrium was 
therefore balanced by pair annihilation. A few synchrotron 
self-compton spectra are shown in Fig.Uand compared with 
previous work. Other output parameters are listed in Table 
[5] Although the general spectrum shape is recovered, sub- 



i - 



B=1000 G 




Fig. 4. Synchrotron self-compton spectra. Spectra are 
shown for 3 magnetic field strengths (Ib = 3.2 x 10~ 2 , 3.2 x 
IP -1 , 3.2), from results of this work (solid lines) and lCoppH 
(|l992j) (dotted lines). Here, l e = 10 and R = 10 14 cm. 



Model 


work 


T e - 




(Ec) 


V3 


1 


NM98 


0.135 


0.085 


0.58 


2.4 




This work 


0.132 


0.082 


0.65 


110 


2 


NM98 


0.080 


0.060 


1.46 


0.44 




This work 


0.086 


0.066 


1.53 


17 



Table 4. Output parameters. r e ± are the electron 
and positron Thomson optical depths. (E c ) = J (7 — 
\)N e -dp/ J N e -dp is the averaged kinetic energy of the 
electron distribution. 7/ = 773 x 10 3 is the heating efficiency 
defined by Eq. [H 



B 


work 


PY 3 




03 


lx/h 


«2-10 


100 


Coppi92 
This work 


21.0 
18.9 


0.508 
0.513 


37.3 
64.1 


82.4 
72.6 


0.578 
0.561 


300 


Coppi92 
This work 


15.4* 
13.5 


0.422 
0.434 


26.9* 
44.7 


58.6 
49.4 


0.609 
0.604 


1000 


Coppi92 
This work 


3.83 
3.29 


0.154 
0.234 


21.4 
28.0 


24.2 
20.4 


0.749 
0.702 



Table 5. Output parameters for synchrotron self-compton 
models. The magnetic field B is given in Gauss. Other pa- 
rameters definitions are the sa me as in Table [U *The pair 
yield and temperature given bv lCoppil (11992ft for B=300G 
were multiplied by 10 to correct for what we believe to be 
typos. 



coefficients found in this work are o ne or two orders of 
magni tude higher than those derived bvlNavakshi n fc Melial 
( 1998j). Given the good agreement in the shape and normali- 
sation of the distributions and escaping spectra, we believe 
that the different efficiencies are the result of a typo in 
their paper. We further ch ecked our e-p Coulomb heating 
rate against the results of iDermer et al.l (|1996t ). and also 
against simple analytic approximations, and found excel- 
lent agreement. 



stantial differences are observed in the far UV, soft X-ray 
band where the flux appears to have been underestimated in 
previous studies. We also derive a measurement of temper- 
ature for the thermal part of the distribution that is larger. 
This most likely results from our more precise treatment of 
the cyclotron emission/absorption. 
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4.4. The synchrotron boiler 

Non-thermal distributions of particles can be thermalized 
by the emission and absorption of synchrotron photons. 
The efficiency of this mechanism however depends on the 
parameters. To i llustrate this p r ocess, we consider the 
case presented in iGhisellini et all (|l998al ). where high en- 
ergy particles strongly emit and absorb synchrotron pho- 
tons. We inject a constant mono-energetic distribution of 
electrons into an empty source of size R = 10 13 cm. 
Electrons escape freely, which leads to a steady state. In 
this study, we consider only cyclo-synchrotron radiation 
and Compton scattering (pair effects, Coulomb scattering, 
and e-p bremsstrahlung are neglected). Figure [5] shows the 
time evolution of the particle and photon distributions and 
a comparison with previous work. 
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Fig. 5. Evolution of the particle distribution (upper panel) 
and the outgoing photon flux (lower panel). Times are 
t/(R/c) = 10- 4 ,10- 3 ,10- 2 ,10- 1 ,1 and 10 from lower to 
higher curves. Particles are injected with a Gaussian distri- 
bution centred at 7 = 10, of width Sj = 1, and with a com- 
pactness parameter: l e - =1. The magnetic compactness is 
Ib = 10 and the domain size is: R = 10 13 cm. The dashed 
curve is the Maxwell-Boltzmann distribution of same nor- 
malization and same average energy as the equilibrium so- 
lution. 



The resul ts qua litatively confirm those of 
Ghiscllin i et al.l (|l998af) . As time evolves, high energy 
particles are cooled by synchrotron radiation, which starts 



to create a radiation field. Soft synchrotron photons 
are up-scattered by high energy electrons and form the 
high energy part of the spectrum. For the choice of 
parameters (l e - /Ib « 1), the effect of synchrotron 
self-absorbed radiation on the particle distribution 
dominates over Compton scattering so that the addi- 
tional cooling on particles by the latter is negligible. 
The synchrotron cooling timescale for particles is then 
t s /(R/c) - ( 7 -l)/(4? B p 2 /3) w l/(7+l)/Jfl and ranges be- 
tween 0.5 Hb for low energy particles and 0.05//b for high 
energy particles (7 ks 20). For Ib — 10, the distribution has 
reached a quasi-thermal shape at f k 0.01 — 0.1i?/c, i.e. on 
the synchrotron timescale (see Fig. [5]). The normalization 
then saturates as the escape of particles balances the 
injection rate, which occurs on a typical time scale of 
t-esc ^ R/c. The low energy part of the distribution is 
well reproduced by a Maxwell distribution but the high 
energy part of the distribution declines more rapidly than 
a real thermal distribution. In more detail, the results 
differ however quiet significantl y. For th e sake of clarity, 
the results of Ghi sellini et al.l fl998a) have not been 
overplotted, athough their distributions are systematically 
colder and broader than the previous ones, especially in 
the transient phase. Although the low energy part still 
looks thermal, the deviation at high energy is therefore 
higher. 

To illustrate the effect of the magnetic field intensity 
and the particle injection rate, we plot in Fig. [5] the tem- 
peratureQ of the steady distribution as a function of the 
injection compactness parameter for 3 different magnetic 
compactness parameters. The steady temperature is quite 
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Fig. 6. Effect ive electron temperatu re estimated following 
Eq. (18) in IGhisellini et all (|l998al ). The 3 curves are for 
magnetic compactness parameters Ib — 0.1, 1,30 from the 
top curve to the bottom curve and the do main size is R = 
10 13 cm . The stars indicate the results of IGhisellini et ail 
(11998aTl for l B = 30. 



insensitive to the magnetic parameter since the tempera- 
ture varies by only a factor of less than 3 as Ib varies by 
over more than 2 orders of magnitude. In contrast, it is 
quite dependent on the injection compactness parameter. 



the temperature is t he effe ctive temperature defined by Eq. 
(18) in IGhisellini et all (|1998al ) 
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For a given magnetic field, the low injection rates produce 
steady states in which the Thomson optical depth is low. 
The Compton cooling is negligible and the final tempera- 
ture is high. As the injection parameter increases, the opti- 
cal depth becomes large, and the Compton cooling becomes 
efficient and dominates at high energy, eventually produc- 
ing far smaller temperatures. The injection of particles at 
an energy far higher than the average particle temperature 
causes the formation of a hard non-thermal tail and a larger 
deviation from the Maxw ell distribution at high energies. 
Compared to the results bv lGhisellini et all (|l998al ). we find 
temperatures that are significantly higher (up to a factor 3) 
at large optical depths. This is probably due to a more pre- 
cise treatment of th e radiat ion field and Compton scatter- 
ing. iGhisellini et alj (|l998a| ) considered only the cooling of 
particles by inverse Compton scattering and assumed that 
it was limited in the Thomson regime. By using the exact 
Klein-Nishina cross section, we find more rapid photon es- 
cape and a weaker radiation field, whose cooling efficiency 
is lower. 

At large Thomson optical depth (i.e. at high injection 
rates), Coulomb exchange is supposed to dominate over 
synchrotron self-absorption. To investigate this, we com- 
pleted the same simulations including e-e Coulomb scatter- 
ing. Results are shown in Fig. \7\ It is found that the e-e 



energy particles becomes very efficient and it is found that 
this effect is significant (up to a factor of 2 for l e = 100). 

A more detailed study of the synchrotron boiler mecha- 
nism and its application to X-ray binaries will be addressed 
in future work. 



4.5. Particle acceleration 

As a second example, we investigate the effect of Fermi, sec- 
ond order acceleration. We consider a magnetised (Zg = 1), 
isolated plasma of size R — 5 x 10 7 cm (typical of X-ray 
binary coronae), with no injection of seed photons. The 
soft photons are emitted by synchrotron radiation of high 
energy particles. The acceleration is modelled by the sec- 
ond order Fermi process and no particle is injected into 
the plasma. Particles are assumed to be trapped and the 
Thomson optical depth is set to be r e = 1. Pair produc- 
tion/annihilation and Coulomb collisions are neglected to 
focus on the role of particle acceleration. After a transient 
phase that depends on the initial conditions, particles and 
photons reach a steady state that depends only on the ac- 
celeration properties. 

We first investigate the role of the acceleration efficiency 
and the threshold energy is assumed to be far lower than 
the bulk of particles. Figure [8] presents the steady particle 
distributions and spectra for various values of the acceler- 
ation efficiency. In all cases, the distribution is similar to a 
Maxwell-Boltzmann distribution. As found in previous cal- 
culations, the diffusion in the m omentum space produces a 
quasi-thermal distribution (e.g. iKatarzvnski et aT] l2006bT ) 
and in this case, the thermalisation is also helped by the 
synchrotron boiler mechanism. The spectrum is the sum 
of the low energy synchrotron emission and a hard tail re- 
sulting from the multiple Compton scattering of these soft 
photons from the highest energy particles. As the accelera- 
tion efficiency increases the steady distribution widens and 
moves to higher energies. As a consequence, the spectra ex- 
hibit a stronger hard tail. The temperature of the distribu- 
tion is given in Table [5] and as expected, it becomes higher 



10 U i i i , i , , , , 

0.1 1-0 10.0 100-0 

Fig. 7. Effect of Coulomb cooling: same as Fig. [6] for Is — 
30, with and without e-e Coulomb scattering (solid and 
dashed line respectively). 



Coulomb collisions tend to increase the effective temper- 
ature. As explained before, particles are injected at high 
energy. They are cooled by both synchrotron emission and 
Compton scattering, and form a low energy thermal pool, 
high energy particles are then scattered by thermal elec- 
trons with e-e Coulomb collisions. The cooling of the high 
energy distribution is very efficient but the thermal pool of 
cool electrons gain energy by this interaction, giving higher 
effective temperatures. This effect is negligible at low in- 
jection rates when the temperature is so high that the in- 
jection energy has a value that is almost in the bulk of the 
distribution and there is no well-marked high energy tail. 
However, at high injection rates, the temperature decreases 
and particles are injected at far higher energies than in the 
thermal pool. Exchange of energy between high and low 



log(/acc) 


tacc 


03 


-4 


48800 


164.5 


-1 


75.1 


259.8 


1 


1.02 


369.0 


2 


0.104 


390.5 



Table 6. Outputs parameters for Fig. [8] in CC is the acceler- 
ation time in unit R/c and 9% is the averaged temperature 
estimated as in section |4~T1 in units of I0~ 3 m e c 2 . 



as the acceleration efficiency increases. Since there are more 
high energy particles, the synchrotron self-absorbed emis- 
sion is stronger, which cools the softer particles more effi- 
ciently and the averaged temperature saturates. The source 
luminosity however scales with the acceleration compact- 
ness parameter. 

We now investigate the role of the minimal energy above 
which the acceleration occurs. The simulations are com- 
pleted with the same parameters but the acceleration effi- 
ciency is set to Z acc = 10 and we vary the threshold energy. 
We note that for a given power supplied to the plasma by 
accelerating particles, the higher the threshold energy the 
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Fig. 8. Particle distributions (upper panel) and spec- 
tra (lower panel) for different acceleration efficiencies. 
Acceleration is modelled by second order Fermi process 
with no threshold. The optical depth is r e = 1, the do- 
main size is R = 5 x 10 7 cm and the magnetic compactness 
parameter is Ib = 1- 




= hf/rr 



Fig. 9. Particle distributions (upper panel) and spectra 
(lower panel) for different threshold energies. Acceleration 
is modelled by second order Fermi process with Z acc = 10. 
The optical depth is r e = 1, the domain size is R = 5 x 10 7 
cm and the magnetic compactness parameter is Ib = I- 
Dashed vertical lines show the threshold momenta. 



fewer the accelerated particles and the shorter their acceler- 
ation time. The results are shown in Fig.[5]for threshold mo- 
menta p c = 0, 1, 1.3, 2, 5, and 10. The corresponding accel- 
eration times are t &cc = 1.02, 6.32 x 10 _1 , 3.45 x 10 _1 , 1.77 x 
10" 1 ,5.63 x 10~ 2 , and 1.95 x lQ- 2 R/c respectively. For 
low threshold energies (p c < 1), the distribution depends 
weakly on the precise threshold. It produces a small num- 
ber of soft photons. Since the particles are not energetic, 
the photons undergo multiple Compton scattering and form 
the strong high energy part of the spectrum. When only 
mid-relativistic or relativistic particles are accelerated, they 
form a high energy tail that extends far beyond the ther- 
mal pool, and the situation then differs significantly. Since 
the total synchrotron emission increases with the energy of 
the emitting particles (/ j s du oc p 2 ), the synchrotron bump 
is much larger. The hard energy tail of the particle distri- 
bution also has a flat slope, which produces a wider syn- 
chrotron bump. Since the particles have higher energy, the 
Compton up-scattering of these soft photons becomes more 
efficient. In the limit where the synchrotron soft photons 
have low energy (w S ynch << 1) and the accelerated elec- 
trons have high energy (7 >> 1), the photon energy gain 
during one single scattering is: w compt /w sync h ~ 47 2 /3. As 
a result, photons undergo only a small number of Compton 



scatterings before they reach the particle energy, forming 
a double-humped spectrum as those of blazars in the case 
p c = 10. At the same time, Compton scattering cools the 
thermal particles further so that the bulk of particles moves 
to lower energies. Future work will consider this effect in 
more details including the physics of wave-particle interac- 
tion. In particular, more consistent models of particle es- 
cape and momentum diffusion must be implemented. 

Conclusion 

We have presented a code developed to model radiation 
processes in high energy plasmas without any assumption 
about the shape of the particle distribution. The code 
is time dependent. It uses the exact Compton and pair 
production/annihilation unpolarized, isotropic cross sec- 
tions. Cyclo-synchrotron self-absorbed radiation is taken 
into account from the sub-relativistic regime to the ultra- 
relativistic one, which represent an improvement in com- 
parison with other codes. It also includes an approximate 
treatment of e-e and e-p Coulomb exchange and e-p self- 
absorbed bremsstrahlung radiation. Explicit prescriptions 
for particle acceleration have also been implemented. The 
code deals consistently with all processes over wide ranges 
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of energy. There is no restriction on the photon energy and 
particles can have momenta in the range 10 -7 < p < 10 7 . 
It can therefore be used to model various sources, such as 
not only X-ray binaries and AGN, but also 7-ray bursts 
and pulsar wind nebulae. 

The major limitation of the code is its simplified ge- 
ometry. The code simulates a uniform system, typically a 
homogeneous sphere with an isotropic and unpolarized ra- 
diation field. It obviously introduces a bias in simulations 
of X-ray binaries coronae, where the seed photons from the 
disc have an isotropic distribution or in jets of AGN, where 
geometrical effects are important. However, the geometry 
of the emitting regions in high energy sources is poorly con- 
strained and in most cases it does not play a crucial role. 
The prescriptions used for particle and photon escape are 
also able to reproduce the main effects of geometry. 

Some examples have been shown of checks of the code 
capabilities in comparison with previous codes designed to 
solve restricted problems, involving a limited number of in- 
gredients. In several cases, we have disabled some processes 
in our code to ensure more rigourous comparisons. We have 
found that the code confirms qualitatively all previous re- 
sults. After considering more precisely these processes, the 
properties of the exact spectra and particle distribution 
were however found to be slightly different. As an example, 
we investigated the acceleration by second order Fermi-like 
processes. We found that an energy threshold for accelera- 
tion produces a non-thermal population of particles when 
it reaches the mid-relativistic regime. 
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